A Fast, High-Order Method in Two Dimensions: Numerical Implementation

In this chapter, we present several significant improvements to the original numerical implementation of the two-dimensional method introduced in [13], as described in the previous chapter. The numerical solution of the associated approximate integral equation (2.3) consists of two main parts: efficient, high-order numerical quadrature rules and an efficient linear solver. The numerical evaluation of the integrals in (2.6) requires both angular and radial integration. In Section 2.3, we described a method for computing the required angular integrals (2.9), exactly by means of Fourier smoothing of the scatterer. Furthermore, this method requires only O(M log M) operations for each radial point. (For smooth scatterers, however, the angular integrals can instead be computed efficiently and with high-order accuracy by direct application of the trapezoidal rule without replacing m with m 2M (see Section 3.1.1).) On the other hand, to compute the radial integrals (K v M)(a) (2.2), we present an improved scheme based on Chebyshev polynomials approximation. More precisely, we approximate I (r) in (2.9) with Chebyshev polynomials on each of several subintervals (see Section 3.1.2). When the integration domain contains the origin, we resolve the logarithmic singularity in the Hankel function by integrating by parts (see Section 3.1.2.1). Then by means of appropriate pre-computed integral moments (see Section 3.1.2.2), we obtain a high-order approximation of the required integrals. Computation of these integrals requires only O(N) operations.